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The equilibrium vortex formations in rotating binary Bose gases with a rotating frequency higher 
than the harmonic trapping frequency is investigated theoretically. We consider the system being 
evaporatively cooled to form condensates and a combined numerical scheme is applied to ensure the 
binary system being in an authentic equilibrium state. To keep the system stable against the large 
centrifugal force of ultrafast rotation, a quartic trapping potential is added to the existing harmonic 
part. Using the Thomas-Fermi approximation, a critical rotating frequency Q c is derived, which 
characterizes the structure with or without a central density hole. Vortex structures are studied 
in details with rotation frequency both above and below Q c , and with respect to the miscible, 
symmetrically-separated, and asymmetrically-separated phases in their non-rotating ground-state 
counterparts. 



I. INTRODUCTION 

Quantum coherence has enabled intriguing phe- 
nomenon, such as quantized vorticity, in Bose-Einstein 
condensates (BECs). When a trapped condensate is 
driven to rotate, singly quantized vortices form. In lower 
rotations, only one or few vortices will be present at equi- 
librium pQ. Faster rotation can generate more vortices 
which are eventually condensed into a lattice [2H5]. Vor- 
ticity in a single-component Bose condensates has indeed 
been observed in a variety of experiments. On the other 
hand, since the first experiment of two coexisting con- 
densates of two different hyperfine states of 87 Rb [6] was 
realized, BEC in mixtures of trapped quantum gases pro- 
vides a unique opportunity to study the miscibility of 
interpenetrating quantum fluid. Several theoretical arti- 
cles about binary-mixture condensates have expounded 
that both the interspecies and intraspecies interactions 
play an important role in determining the density pat- 
terns and phase separation of the condensates j7MT2] . In 
contrast to the Abrikosov vortex-lattice state of a scalar 
BEC, the vortex states of binary- mixture BECs have var- 
ious exotic structures due to the variety of interactions 
H3HT5I. 



Rotating Bose condensates are usually confined in a 
harmonic trap with cylindrical symmetry around the ro- 
tation axis (say, z-direction). In these typical cases, there 
are two limiting regimes depending on the relative size of 
the rotating frequency and the trapping frequency uj 
in the xy plane. When > uj, the system will become 
unstable due to a strong centrifugal force. In order to an- 
alyze the regime of ultrafast rotations with > uj, one 
approach is to add a quartic part to the harmonic poten- 
tial. In this type of system, the trapping force will be al- 
ways greater than the centrifugal force and consequently 
the regime > uj can be fully explored [T6H23] . The 
current paper attempts to study the equilibrium vortex 
states of ultrafast-rotating binary condensates confined 
in a harmonic-plus-quartic potential. 

In a single-component fast-rotating condensate, when 
Qo > uj, the system can experience an effective potential 



of Mexican-hat shape. Depending on how Qq is larger 
than cj, the system can be roughly separated into two 
regimes: a condensate with or without the central density 
hole. More exactly, it has been shown in the literature 
that there exist three distinct phases for the fast-rotating 
scalar condensate confined in a harmonic-plus-quartic po- 
tential. One is the vortex lattice without a hole (VL), the 
second is the vortex lattice with a hole (VLH), and the 
third is the giant vortex state (GV) [19j [22] . It is inter- 
esting to see how the interspecies and intraspecies inter- 
actions play the role in the binary-mixture condensates 
under fast rotation and in particular, how the above- 
mentioned phases manifest in these systems. 

Due to the complexity of the interactions in the bi- 
nary system, a standard imaginary-time propagating 
method may not be easy to find the converging re- 
sults for the true equilibrium states. It has been shown 
in Ref. [24] that Gross-Pitaevskii equation (GPE) with 
a phenomenological damping term can provide an effi- 
cient numerical machinery for finding the eigenstates of 
the time-independent GPE. This, a similar stochastic 
Gross-Pitaevskii equation (SGPE) approach [25ti29] , has 
been demonstrated to be an efficient way for studying 
the single-component BEC system. It is also anticipated 
that SGPE is an alternatively efficient method for study- 
ing the dynamic and equilibrium properties of a multi- 
component system near absolute zero. 

The paper is organized as follows. In Sec. [TTJ we in- 
troduce the theory for studying the equilibrium vortex 
states of a binary-mixture BEC system. To investigate 
the regime of ultrafast rotation, the system is trapped 
under a harmonic-plus-quartic potential. Sec. |III| is de- 
voted to a detailed discussion of the vortex structures 
for fast-rotating binary-mixture BECs. A critical ro- 
tation frequency Q c is derived and both Q < Q c and 
Q > Q c regimes are studied. It will be shown explic- 
itly that vortex structures of the system do manifest the 
ground states of their non-rotating counterparts. Sec. [TV 
is a brief conclusion. 
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II. METHODOLOGY 

We consider rapidly rotating two-component pancake- 
shape BECs that are parallel to the x?/-plane and in a 
cylindrically symmetric potential. Assuming that the ex- 
citation in the z-direction is suppressed, the system can 
be treated approximately by a two-dimensional theory. 
In the mean-field approximation, such a two-component 
BEC system in a co-rotating frame with a rotating fre- 
quency around the z-axis can be described by the 
time-dependent coupled GPEs: 
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where for component j (= 1, 2), Cqp is the GP Hamilto- 
nian, \i- is the chemical potential, and ^ 3 is the macro- 
scopic wave function normalized under Nj = J \^/j\ 2 dxdy 
with Nj the particle number. Hj = —h 2 V 2 / (2m j) + 
Vj — ^0^2 is the single-particle Hamiltonian with rrij 
the atomic mass, Vj = rrijUJ 2 r 2 /2 + Ujr A /A the trap- 
ping potential in polar coordinates (r, 0), cjj the har- 
monic trapping frequency, the strength of the quar- 
tic potential, and L z = —ihd/dcj) the z-component an- 
gular momentum operator. The interaction parameter 
Ujk = 2irh 2 ajk ( m ^j l + m k 1 ) w ^h djk(> 0) the effective 
two-dimensional s- wave scattering length between atoms 
in components j and k. 

As mentioned before, it has been shown that GPE with 
a phenomenological damping term [24 and SGPE may be 
an efficient numerical approach to obtain accurate ground 
states of a given time-independent GPE. Here we shall 
apply a similar SGPE approach to study the vortex states 
of a ultrafast-rotating two-component BEC system. The 
coupled SGPEs for the present binary system can be ex- 
pressed as 



(2) 



where r} 1 = 77 (r,t) is a complex Gaussian noise 
considered arising due to the contact with the ther- 
mal modes. The correlation function associated 
with the noise is given by (77* (r, t) r} 1 (V, = 
2hk B TjS (r - r') S (t - t'). The strength of the noise, 
and hence the damping, is thus proportional to 7. In 
principle, 7 can be calculated ab initio in terms of 
the Keldysh self-energy [25j [26]. However, as we are 
interested only in the properties at equilibrium, we 
may approximate it as a spatial and temporal constant. 
Throughout this paper, 7 is reasonably taken to be a 
small number, 0.01. The temperature is set to be InK 
which is about 10 -2 or less of the critical temperature 
of a typical BEC system. Here we have to state that 
the coupled SGPEs ([2| do not describe the dynamics of 
a two-component system accurately since any coupling 
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FIG. 1. Phase diagram of the non-rotating binary-mixture 
condensates confined in a harmonic-plus-quartic trap in terms 
of the relative interaction strengths, oli = 022/011 and ai2 = 
0*12 /o>ii- The quartic trap strength is fixed at A = 1. Re- 
gions I, II, and III correspond to miscible, asymmetrically 
separated, and symmetrically separated phases, respectively. 



of the thermal cloud components of the two different 
condensates is not considered. Because we focus on the 
equilibrium states, the effect of these couplings should 
be slight, and Eq. Q is good enough for obtaining the 
equilibrium states of a two-component system. 

To reduce the number of parameters, we shall assume 
that mi = rri2 = m, uj\ = U2 = cj, u\ = U2 = u, 
and Ni = AT 2 = iV, respectively. Moreover, for conve- 
nience, the computations will be carried in the oscillator 
units. That is, the lengt h, tim e, and energy are scaled 
respectively in units of y/h/muj, 1/cj, and fouo. As a con- 
sequence, the coupled SGPEs ([2| take the following di- 
mensionless forms: 
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Here we have redefined the normalized wave function 
ipj = y/h/ (muoN)^ j, the strength of the quartic trap 
A = uhj (m 2 cj 3 ), the interaction constants between 
atoms gjk = AirNdjk, the chemical potential jlj = /ij/huj, 
and the noise 77 = tj^/Huj. Besides the rotation rate 
Q = Qq/lj. The rotation rate and the quartic trap 
strength will be fixed at Q = 2.5 and A = 1 in our calcu- 
lation throughout this paper. 

In our calculations, the equilibrium solutions and the 
corresponding chemical potentials are obtained by solv- 
ing the norm-preserving imaginary-time propagation of 
the time-dependent coupled GPEs ([!]) starting from an 
arbitrary trial wave function. The propagation contin- 
ues until the fluctuation in the norm of the wave func- 
tion becomes smaller than 10 -5 . To determine whether 
the vortex states obtained are indeed at equilibrium, the 
solutions of the GPEs, which were converged by the 



imaginary- time method, are then substituted into and 
treated as the initial states of the coupled SGPEs If 
the initial state was not an equilibrium state, it would 
keep propagating until the damping term vanishing. 

Moreover, we have used the method of lines with spa- 
tial discretization by the Fourier pseudospectral method 
to compute Eqs. ([!]) and Q [or The time integration 
in Eq. ([!]) is done by the adaptive Runge-Kutta method 
of order 2 and 3 (RK23), which is more time efficient due 
to an adjustable time step. However, the fourth-order 
Runge-Kutta method (RK4) is used for Eq. Q [or p^, 
since the thermal noise term r] 1 is proportional to 1/ Vat, 
and is better computed with a fixed time step. 



III. RESULTS AND DISCUSSIONS 

A. Non- rotating ground states 

Solving the non-rotating (Q = 0) time-dependent cou- 
pled GPEs using the imaginary-time propagating 
method, we have obtained three distinct phases for the 
ground states of the binary-mixture condensates. Fig. [I] 
shows the phase diagram of the binary-mixture con- 
densates in terms of the relative interaction strengths, 
a 2 = ^22/^11 and a 12 = By symmetry, it is 

sufficient to consider 0L2 < 1 only. The non-rotating 
binary-mixture condensates are phase miscible in region 

I and phase separated in regions II (asymmetric) and 
III (symmetric). As shown in Fig. [I] the boundary be- 
tween phases I and III is linear, ql\2 = 0^2-, which is ob- 
tained by jointing nine points: from a 2 = 0.2 to ct2 = 1 
spaced by 0.1, and theirs corresponding ol\i are deter- 
mined numerically with precision less than 0.01. The 
quartic trap (A) has only a minimal effect on the bound- 
ary between the miscible and the symmetric separated 
phases. On the contrary, the boundary between phases 

II and III is quite A-dependent, which is obtained by 
jointing nine points: from a 2 = 0.4 to ct2 = 1 spaced by 
0.1, and theirs corresponding ol\i are determined numer- 
ically with precision less than 0.01. Due to a relatively 
stronger confinement, the area of phase II can expand 
as A increases. This allows a larger space for studying 
the asymmetric phase-separated regime for fast-rotating 
binary-mixture condensates. The two boundaries inter- 
sect at (0^2,^12) = (1? 1) which is called the isotropic 
point. We shall study the vortex structures in all three 
phases and in particular at the isotropic point. 



B. Vortex states 

Before discussing the fast-rotating binary-mixture con- 
densates, it is useful to first examine the case of a scalar 
rotating condensate trapped in a harmonic-plus-quartic 
potential. In the one-component system, there is only 
one interaction constant g = 47riVa, and when g ^> 1 the 
system is in the so-called Thomas-Fermi (TF) regime. 
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FIG. 2. (Color online) Vortex structures of individual compo- 
nent in fast-rotating binary-mixture condensates confined in 
a harmonic-plus-quartic potential. The interaction constants 
9n = 922 — 1300 {a.2 — 1), the rotation rate ft = 2.5, and the 
quartic trap strength A = 1 for all frames. From the top to 
the bottom rows, olyi — 0.7, 1, and 1.3, respectively. 



The TF density can be obtained by ignoring the kinetic 
energy: 



n TF |^tfI 



/i TF 



fo 2 - 1) r 2 Ar 4 
2 4 



/g- (4) 



For Q < 1, the density has a local maxium near the 
center, but it changes to a local minium for Q > 1. Eq. 
Q has solutions where the TF density vanishes 
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(5) 



Here the upper (plus) sign denotes the outer radius 
R > for any value of the chemical potential /i XF . In 
contrast, the lower (minus) sign yields a physical in- 
ner radius R < only if /i TF is negative. Correspond- 
ingly the negative chemical potential can be fixed by 
/i TF = (3^v^A/tt) 2 / 3 /4 - (ft 2 - 1) 2 /(4A) and the spa- 
tial extension of the atomic cloud can be characterized 
by R> and R < which correspond to the outer and in- 
ner TF radii, respectively. Accordingly, the system can 
be separated into two regimes demarcated by a critical 
rotating frequency Q c corresponding to R < (Q C ) = [or 
/i TF (fi c ) = 0]. It is found that ft 2 = 1 + (3X 2 g/7r) 1 ^ 
which is dependent of the interaction constant g and the 
quartic trap strength A. At sufficiently high rotating fre- 
quencies such that Q > Q c and hence > 0, a central 
hole will appear in the condensates, i.e., in the VLH state 
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FIG. 3. (Color online) Vortex structures of component 1 (left 
panel) and 2 (right panel) in fast-rotating binary-mixture con- 
densates. Here gu = 1300, ct2 = 0.5, olyi — 0.8, A = 1, and 

n = 2.5. 



P~9j|22]. In contrast for Q < Q c , vortex lattice without a 
central hole regime (the VL state) will appear. 

For the present binary-mixture condensates, the criti- 
cal rotating frequency Q c which characterizes the transi- 
tion between the VL and the VLH states, can be qualita- 
tively determined by the larger of gu and g 22 . As men- 
tioned before, in this paper we consider only the cases 
9n > 922 (i-e., a 2 < !)• Consequently 



EE 1 + (3A 2 0n/7r) 



(6) 



When Q < Q c , the rotating binary system will be in the 
VL state, while when Q > Q c , the system will be in the 
VLH state. In the following two subsections, two dis- 
tinct cases of g\\ = 1300 and gu = 55 will be studied. 
The former corresponds to a critical rotating frequency 
Q c = 3.43 and hence Q = 2.5 < Sl c , while the latter 
corresponds to Q c = 2.18 and hence Q > Q c . In the 
following calculations, the equilibrium vortex solutions 
and the corresponding chemical potentials are obtained 
by solving the norm-preserving imaginary-time propa- 
gation of the time-dependent coupled GPEs start- 
ing from an arbitrary trial wave function. To determine 
whether the vortex states obtained are indeed at equilib- 
rium, the solutions of the GPEs, which were converged 
by the imaginary-time method, are then substituted into 
and treated as the initial states of the coupled SGPEs 



1. VL state with Q < Q c 

We first consider the vortex structures with rotation 
frequency below the critical rotating frequency, ft < 
Q c . Fig. [2] shows the equilibrium vortex structures of 
two-component condensates confined in a harmonic-plus- 
quartic potential with gu = g 22 = 1300 (a 2 = 1) and 
ol\2 = 0.5, 1, and 1.3 (from the top to the bottom) re- 
spectively. All three cases are belonging to the VL regime 
to which one is able to conclude for the following, (i) 
For a phase-miscible mixture, the equilibrium state is 
composed of regular vortex lattices which form roughly a 
square lattice at a± 2 = 0.5 than what is expected to be a 
triangular lattice when a\ 2 — >• [30]. (ii) At the isotropic 



point (a 2 = ol\ 2 = 1), a honeycomb lattice is formed 
for one component, while vortices in the other compo- 
nent form a vortex-pair lattice (vortex of every pair has 
the same circulation), (iii) Stationary vortex sheets are 
formed for an asymmetric phase- separated mixture. Our 
results in Fig. [2] are intended to be compared to those 
shown in Figs. 2(a), 3, and 4 in Ref. [14]. With our re- 
sults, we have been able to verify that the one shown in 
Fig. 3(b) of Ref. [14] corresponds to an authentic equilib- 
rium state, while the one shown in Fig. 3(a) of Ref. [T4] 
corresponds to a transition state. 

Fig. [3] shows the vortex structures of the mixtures with 
#11 = 1300, a 2 = 0.5, and a\ 2 = 0.8. For these param- 
eters, the system is still in the VL regime with Q < O c . 
In this case, the system has a ball-and-shell non-rotating 
ground state (phase III) and in the vortex state it forms 
an interlocking oniony vortex-sheet structure. 



2. VLH state with Q > Q c 

The results in the VLH regime with Q > Q c are con- 
sidered next. In this subsection, to see more clearly the 
vortex physics, we show both results of both density pro- 
file rij(x,y) = \ifjJx,y)\ 2 and phase profile given by 



Sj(x, y) = arctan 



"km/)Ax,y) 



ReipJx.y) 



(7) 



In the phase profile, the end point of the boundary be- 
tween a 7r phase line and a — ir phase line will correspond 
to a vortex. In addition, the circulation and the number 
of vortices can also be counted directly. Fig. [4] shows the 
vortex structure and the corresponding phase profile of 
fast-rotating binary- mixture condensates with gu = 55, 



4. 
21 

^ ol 

-2 1 



O 



0.08 
0.06 
0.04 
1 0.02 



4. 
2 1 

^ ol 

-2 1 



O 



0.08 
0.06 
0.04 
0.02 




FIG. 4. (Color online) Vortex structures (top row) and phase 
profiles (bottom row) of component 1 (left column) and 2 
(right column) of fast-rotating binary-mixture condensates. 
Here gu — 55, ct2 — 1, ai2 — 0.5, A = 1, and Q — 2.5. 
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FIG. 5. (Color online) Vortex structures (top row) and phase 
profiles (bottom row) of component 1 (left column) and 2 
(right column) of fast-rotating binary-mixture condensates. 
Here gu — 55, a 2 — 1, olyi — 1.3, A = 1, and Q — 2.5. 



^2 = 1, and ai2 = 0.5. The ground state of the cor- 
responding non-rotating condensate mixture is miscible 
(phase I) to which the wavefunctions of the two compo- 
nents overlap entirely. In view of Fig. [4| it is found that 
the two annular vortex arrays are interlocking in a man- 
ner that density peak of one component is located at the 
density hole of the other component. 

Fig. [5] shows the vortex structure and the correspond- 
ing phase profile of the mixture with gu = 55, a 2 — 1, 
and ol\2 = 1-3. The non-rotating counterpart has an 
asymmetric separated ground state (phase II in Fig. [I]). 
Due to the strong repulsive interaction between two con- 
densates which results the asymmetric separated charac- 
teristic, the two condensates occupy on the opposite side 
of each other. 

Fig. [6] shows the vortex structure and the correspond- 
ing phase profile of a condensate mixture with gu = 55, 
OL2 = 0.5, a\2 = 0.8. With respect to phase III in 
Fig. [I] the non-rotating mixture has a ball-and-shell 
ground state, i.e., the component with larger intraspecies 
interaction occupying the outside and forming a shell, 
while the component with smaller intraspecies interac- 
tion occupying the inside and forming a ball. It is found 
that vortices in component 1 (of larger intraspecies in- 
teraction) form a circular array around the central low- 
density hole, while vortices in component 2 (of smaller 
intraspecies interaction) also form a circular array which 
is closer to the center (see the inset of the density profile 
or the phase profile plot). In addition to the vortices, 
the "ball" of the non-rotating counterpart of component 
2 is actually pushed away (due to fast rotation) from the 
center and forms a robust ringlike condensate located at 
where the vortices of component 1 are (i.e., interlocking). 
This robustness actually resists its own vortices revolving 
into it. 

In Fig. [7| we have also shown the vortex structure 




FIG. 6. (Color online) Vortex structures (top row) and phase 
profiles (bottom row) of component 1 (left column) and 2 
(right column) of fast-rotating binary-mixture condensates. 
Here gu 55, a 2 0.5, a 12 0.8, A = 1, and Q = 2.5. The 
inset shows a more clear view of the vortices in component 2. 



and its corresponding phase profile of a mixture at the 
isotropic point, with gu = 55, oli = <^i2 = 1- Similar to 
the result of the middle row of Fig. [2j fast -rotating con- 
densates tend to form the vortex-pair structure at the 
isotropic point. 



C. Concluding Remarks 

The overall features of vortex structures in ultrafast- 
rotating binary-mixture condensates can be understood 
as follows. Pertaining to phase I in Fig. [I] when a\2 = 
where the two components are not interacting with each 
other, the theory is essentially reduced to the one for 
single components. In this limit, triangular vortex lat- 
tices are expected to form with Q < O c , while annular 
vortex arrays are expected to form with Q > Q c . As 
a\2 is present and increases, vortex cores of one compo- 
nent gradually shift away from those of the other com- 
ponent and consequently with Q < Q c the triangular 
lattices are distorted. Eventually the vortices for each 
component will form a square lattice instead of a trian- 
gular one. As a 12 exceeds «2, equivalently for the sys- 
tem to shift to phase II or III, the condensates can un- 
dergo phase separation to spontaneously form domains. 
For one condensate, the cavity of another condensate 
is where the lower effective potential is, which is more 
apt to be occupied. However, more interlocking will 
cause wavefunctions overlap more and at the same time 
raise the interspecies interaction energy oc gi2 |^i| 2 1^2 1 2 - 
In order to prevent the above-mentioned interlocking 
that causes high energies, the vortices are actually con- 
centrated out of the condensates and form the vortex 
sheets. Consequently the two components form comple- 
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FIG. 7. (Color online) Vortex structures (top row) and phase 
profiles (bottom row) of component 1 (left column) and 2 
(right column) of fast-rotating binary-mixture condensates. 
Here gn — 55, ol^ — 1, olyi = 1, A = 1, and Q — 2.5. 



trapped in a harmonic-plus-quartic potential. In contrast 
to the harmonic trap alone case where the system is un- 
stable when the rotation frequency is higher than the 
radial trap oscillator frequency cj, the added quartic trap 
can lead the system to remain stable at higher rotation 
velocity (Oq > lj). Due to the complexity of interactions 
in the binary system, there often occur many metastable 
states in the fast-rotating two-component condensate sys- 
tem and the standard imaginary-time propagating ap- 
proach may not really converge to the true equilibrium 
states of the system. In this regard, we have applied a 
combined numerical scheme to effectively assure that the 
density profiles do really saturate at sufficiently low tem- 
peratures. A critical rotating frequency Q c which charac- 
terizes the transition between the VL and the VLH states 
is identified. Under high rotation frequencies (Qq > cj), 
a variety of vortex structures of the two-component con- 
densates are shown for Q < Q c , similar to those presented 
in Ref. [H] , and also for Q > Q c in particular to which 
various annular vortex structures occur. 



mentary structures to each other and the total density 
will be roughly described by the Thomas-Fermi distribu- 
tion |^i| 2 + |^ 2 | 2 oc max[/i TF - (muj 2 r 2 /2 + ur 4 /4),0]. 
The results in the bottom row of Fig. [2] as well as in 
Figs. 3, 5, and 6 are examples of this kind. 

IV. CONCLUSION 

This paper investigates the equilibrium vortex struc- 
tures of ultrafast-rotating binary-mixture condensates 
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